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Abstract 

A method for computing Ice shapes along the 
leading edge of a wing and a method for predicting 
Its aerodynamic performance degradation due to 
Icing Is described. Ice shapes are computed using 
an extension of the LEWICE code which was developed 
for airfoils. The aerodynamic properties of the 
Iced wing are determined with an Interactive scheme 
In which the solutions of the Invlscld flow equa- 
tions are obtained from a panel method ‘and the 
solutions of the viscous flow equations are 
obtained from an Inverse three-dimensional finite- 
difference boundary-layer method. A new Interac- 
tion law Is used to couple the Invlscld and viscous 
flow solutions. 

The application of the LEWICE wing code to the 
calculation of Ice shapes on a MS-317 swept wing 
show good agreement with measurements. The Inter- 
active boundary-layer method Is applied to a 
tapered Iced wing In order to study the effect of 
Icing on the aerodynamic properties of the wing at 
several angles of attack. 

1 .0 Introduction 

The National Aircraft Icing Technology Plan 1 
called for the creation of a method to simulate ice 
accretion on an aircraft and predict the effect of 
Ice accretion on aircraft performance. This paper 
presents results which suggest that the realization 
of this goal Is within reach, largely due to the 
technical and organizational abilities of NASA 
Lewis Research Center, which developed the LEWICE 
code for the calculation of Ice-accretion on air- 
foils as a function of atmospheric conditions 2 * 3 
and supported the development of methods to deter- 
mine the resulting aerodynamic properties of Iced 
airfoils and wings. At the same time, the 
NASA Lewis Research tunnel has been used for the 
study of the process of Ice formation, and an 
experimental program has been Initiated to measure- 
degradation of the aerodynamic performance of air- 
foils' and wings. 8 It Is Intended that the 
resulting simulation method, which will solve the 
conservation equations and predict Ice accretion, 
will be evaluated by comparing the results of these 
calculations with experimental data. This method 
will then be used to assist In the prediction of 
the performance degradation of wings In Icing con- 
ditions and thereby assist In ensuring that the 
Regulations of the Federal Aviation Authority 9 
are met. 

Research on aircraft Icing can be judged In 
many ways; for example, by the presentations made 
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to the Icing Workshop held at NASA Lewis Research 
Center annually and by the several papers presented 
each year at the AIAA Aerospace Sciences Meeting. 
The experimental Investigation of Bragg et al.' 8 . 
In which the effects of simulated accretions of 
glaze Ice were measured on a wing with a NACA 0012 
airfoil cross-section, follows from a series of 
earlier valuable contributions .”~ 13 Numerical 
investigations have already begun to consider 
fixed and rotary wing performance (l.e. Kwon and 
Sankar 14 * * 8 ) . A major requirement, however. Is 
to extend the LEWICE code to the representation of 
Ice shapes on the leading edge of a wing and 
combine It with a method for the determination of 
performance, and to do so for three-dimensional 
flows over complete airframes and at a modest cost 
so that the Implications of aircraft geometry, 
angle of attack, speed, and atmospheric conditions 
can be predicted. It Is self-evident that this 
achievement will be of value only If the components 
of the simulation method as well as the full simu- 
lation method In Its entirety are tested and the 
results compared with experimental data. 

It Is generally accepted that the variation of 
the shapes of Ice with atmospheric conditions which 
lead to rime Ice can be determined on airfoils with 
the current version of the LEWICE code and that the 
shape of glaze Ice, although acceptable, can and 
should be Improved. 18 *” In this work It Is 

assumed that the LEWICE code Is the best available, 
generally satisfactory, and amenable to Improve- 
ment. Thus, the problem Is to develop a three- 
dimensional version of the LEWICE code and combine 
It with a three-dimensional method for the deter- 
mination of aerodynamic properties, taking Into 
account the fact that Ice formations are neither 
continuous nor smooth. Since one of the main com- 
ponents of the present computer method to predict 
Ice shapes Is a method for the solution of 
Invlscld-f low equations, it would be advantageous 
If this approach could also be used for the predic- 
tion of aerodynamic properties of Iced wings. 

It can be envisaged, at least by those with a 
long-term view, that the Reynolds-averaged Navler- 
Stokes equations might provide the basis for the 
prediction of aerodynamic properties, and scrutin- 
izing the related literature on the solution of 
two-dimensional equations for smooth airfoils 
reveals that there Is a dichotomy of view. It 
would seem, however, that even these comparatively 
simple flows are Investigated using equations that 
are further reductions of the Navler-Stokes equa- 
tions, Including thin-layer Navler-Stokes equations 
and by the Interaction of Invlscld and boundary- 
layer equations. The former approach, like those 
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with higher-order equations, holds much promise, 
but the cost of computation - even In two dimen- 
sions - Is large and becomes prohibitive for use 
In three dimensions, particularly with complete 
aircraft geometries. For the foreseeable future 
the interactive approach Is to be preferred. The 
Invlscld-f low procedure may be Identical to that 
already utilized within the LEWI CE code so that 
the Ice accretion may be part of an Interactive 
calculation. 

The Interactive method used In this paper has 
been widely used for two-dimensional flows, so Its 
prediction abilities are known for a wide range of 
smooth body shapes, Iced airfoils, and aerodynamic 
and atmospheric conditions. In addition, It has 
recently been applied to wings 18 and, although 
further developments are desirable to determine a 
more efficient means of achieving the Invlscld- 
vlscous Interaction In three-dimensional fpows, It 
is clearly able to deal with wing geometries and 
potentially complete aircraft. In this paper, we 
describe a general procedure for predicting Ice 
shapes on wings and determining their Influence on 
the aerodynamic properties of the wings. This 
procedure Incorporates a general 1 nvl sc Id code 
developed by Hess. 19 

This paper describes the extension to wing 
flows of the combined LEWICE/IBL procedure devel- 
oped for airfoils. Section 3 describes the pro- 
cedure developed for calculating Ice shapes along 
the leading edge of a wing. Section 4 describes 
the Interactive boundary-layer method for computing 
three-dimensional flows on Iced wings. In both 
sections, use Is made of Hess's three-dimensional 
panel method for computing the 1 nvl sc 1 d flowfleld 
as described In Section 2. Section 5 presents a 
comparison between calculated and measured Iced 
shapes obtained under Infinite swept wing condi- 
tions on a MS-317 swept wing. The results of Sec- 
tion 6 show the Influence of 390 and 1164-second 
rime Ice on the performance degradation of a NACA 
0012 tapered wing at several angles of attack. In 
this case, the Interactive boundary-layer calcula- 
tions were performed for a prescribed Ice shape 
distribution along the wing leading edge using the 
procedure described In Section 4 . The paper ends 
with a brief discussion of the Implications of the 
work and a summary of the more Important con- 
clusions. 

2.0 Three-Dimensional Panel Method 

He use the surface-source panel method devel- 
oped by Hess 19 for computing the Invlscld flow 
about a complete airplane configuration. This 
method represents a general body by a set of quad- 
rilateral panels as shown In Fig. 1. A three- 
dimensional configuration In general consists of 
lifting sections, such as a wing or pylon for which 
there Is a well-defined trailing edge, and nonlift- 
ing sections, such as the fuselage. Under the 
formulation adopted by Hess, all panels are 
assigned an Independent source distribution, while 
those on a lifting section are assumed to carry a 
bound vortlclty distribution. The variation of 
this bound vortlclty In the streamwlse direction 
Is assumed, while Its variation In the spanwlse 
direction Is adjusted to satisfy the Kutta condi- 
tion at the trailing edge. The complete solution 
for a prescribed flow condition is obtained by 
simultaneously satisfying a condition of zero nor- 
mal velocity at a control point on each panel of 



Fig. 1. Wing/body configuration. 

the body together with a Kutta condition at each 
trail Ing-edge panel . 


The nature of the Kutta condition varies 
greatly among the many panel methods which are 
currently available. While the Kutta condition 
adopted by Hess assumed equal upper and lower sur- 
face pressures at the trailing edge, other methods 
make use of other derived conditions that do not 
guarantee a pressure match at the trailing edge. 
For Instance, Margason et al. 20 showed that 
pressure mismatches of up to half of the free- 
stream dynamic pressure could occur from such 
alternate forms of the Kutta condition. Since we 
are also Interested In the computation of Invlscld 
flows with viscous effects, and since the behavior 
of the boundary-layer at the trailing edge can have 
a significant effect on the overall solution, It 
Is believed that the approach adopted here Is more 
realistic and blends In well with the previous 
studies on two-dimensional flows employing Inter- 
action procedures. 


3.0 Three-Dimensional Ice Shapes and 
Their Calculation 

Before we describe our procedure for calcu- 
lating the Ice shapes along the leading edge of a 
wing for given atmospheric conditions, It Is use- 
ful to review the Ice accretion procedure used In 
the LEWICE program for airfoils. LEWICE consists 
of three main modules: (1) prediction of Invlscld 
flow by a panel method, (2) calculation of the 
trajectories of water droplets driven by the air 
velocities so that the impingement pattern of the 
water droplets on the surface can be determined, 
and (3) a quasi-steady-state surface heat transfer 
analysis In which mass and energy equations are 
solved for a two-dimensional flow In order to 
determine the Ice shape. As discussed In Ref. 3, 
the heat balance equation for a control volume 
Includes terms providing heat to or removing heat 
from the control volume. The dominant "heat 
removal* term Is due to convective cooling and the 
dominant heat source term Is due to the heat 
released by the freezing of water entering the 
control volume. The water entering a control vol- 
ume consists of: (1) water droplets Impinging on 
the airfoil surface, and (2) water *runnlng back" 
from an adjacent upstream control volume. 
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The above procedure for airfoils has been found 
to be remarkably good for predicting Ice shapes, 
(see for example Ref. 17) despite the very 

empirical nature of the expressions used In the 
heat balance as well as In the physical model of 
the Ice accretion process. Studies conducted by 
Olsen and Walker^ *^2 show that at aircraft 

speeds there Is no flow of liquid over the surface 
of the Ice after a short Initial flow, even at 
subfreezing temperatures that are close to the 
freezing point. Instead, there are very large 
stationary drops on the Ice surface that lose 

water from their bottoms by freezing as water Is 

added to the droplets by the addition of 

microscopic cloud droplets. This observation 
disagrees with the existing physical model used In 
the LEWICE code, since It assumes that there Is a 
thin liquid continuously flowing over the ice 
surface. While this model disagrees with the 
experimental evidence, especially for conditions 
leading to horn-shaped clear (glaze) Ice, ‘It does 
reasonably well for rime Ice which forms when the 
air temperature Is low enough to cause the cloud 
droplets to freeze almost Immediately on Impact. 


The role that the liquid layer plays In the 
formation of glaze Ice has been further elaborated 
on by Hansman and Turnock^ and Bllanln^*, who 
point to the surface tension phenomenon as the most 
likely source for the observed deviations. There 
Is convincing evidence from tests with variable 
surface tension that Ice growth depends on TT“ 
Although a smooth water film may form at the stag- 
nation point, It soon breaks up by the beading 
process and freezes Into Irregular rough lumps that 
enhance heat transfer and speed up the Ice growth 
The parameters controlling transition from a smooth 
film to beads and to rough Ice formation are not 
well known at this time. For this reason, care Is 
required to extend the Ice accretion procedure from 
two- to three- dimensional flows. 


In principle, the extension of the LEWICE air- 
foil procedure to three-dimensional flows requires 
that the trajectory calculations be performed for 
the three components of the velocity obtained from 
a three-dimensional panel method and a heat trans- 
fer balance which also Includes the spanwlse 
direction. While the former can be done with 
increased computer cost, the latter Is difficult 
In view of the empirical nature of the heat trans- 
fer coefficient, hr, used to represent the con- 
vective cooling. In the present LEWICE program, 
this parameter for turbulent flows is expressed by 


h c - 1 


c f /2 


Pr t ♦ /(c f /2) (1/St k ) 


] P u e c p (1) 


through the Reynolds analogy with local skln- 
frlctlon coefficient, cr, given by an expression 
valid for a flat-plate flow with roughness, 


If . OM , 

2 * l ln [864 ( 6 t /k s ) ♦ 2.568]’ 


(2) 


The parameter St|< represents a dimensionless rough- 
ness parameter given by 

u k, - 0.2 

St fc - 1.16 (~) (3) 

The momentum thickness, e*. Is calculated by a 
momentum Integral method which employs a smooth 
flat-plate skin-friction coefficient as one of the 


auxiliary equations necessary to solve the Integral 
equation. Since Eq. (1) Is based mainly on data 
correlation with roughness parameters chosen so 
that the computed Ice shapes fit the experimental 
data, any attempt to add a spanwlse heat transfer 
coefficient to the heat balance equation will be 
empirical and can only be done where experimental 
data exists. As a first step, however, It Is pru- 
dent to leave the heat balance In Its two- 
dimensional form and assume It to apply to a three- 
dimensional body which Is expressed In an equiv- 
alent two-dimensional form. One approach, followed 
by Potapczuk and B1dwell“, is to perform the 
trajectory calculations for a three-dimensional 
flowfleld by determining the Impingement patterns 
of the water droplets resulting from the three 
components of the velocity field and apply this 
along the streamlines on the wing. Another 
approach, followed here, Is to approximate the 30 
wing by an equivalent yawed Infinite wing at each 
spanwlse station. In this case, the flowfleld can 
be calculated by a three-dimensional panel method 
and the particle trajectories calculated In a 
quasl-three-dlmenslonal manner, so that the two- 
dimensional heat balance can be applied along the 
airfoil section normal to the leading edge of the 
wing. Each approach has Its advantages and must 
be evaluated against experimental data. Studies 
reported by Potapczuk and Bldwell show that their 
approach leads to Ice shapes which are In good 
agreement with experimental data. 

It Is also useful to seek other approaches to 
extend the two-dimensional Ice accretion model to 
three-dimensional flows since the computer cost of 
predicting an Ice shape at a given spanwlse section 
Is relatively high. For example, the approach of 
Potapczuk takes about 20 minutes on a Cray, while 
the present approach takes less than 7 minutes. If 
the wing is approximated, say by ten airfoil sec- 
tions, then the use of the computer program to pre- 
dict the Ice shapes along the leading edge of a 
wing can require between one and three hours. It 
Is also possible that the dominant effect lies with 
the heat transfer balance and that "more accurate" 
trajectory calculations may not have an appreciable 
effect on the Ice shapes. The answers to all these 
questions, however, can be provided by comparing 
the calculated Ice shapes with experimental data. 


3.1 Present Approach for Calculating Ice Shapes 

on Mings 


To describe our approach, which Is based on the 
extension of the LEWICE airfoil code to wings, we 
use the three-dimensional panel method described 
In Section 2. For a wing defined by streamwlse 
cross-sections, we calculate the velocity compon- 
ents VJ, VJ, V; In a Cartesian coordinate system x, 
y, z at a specified Incidence angle a measured from 
the x-axIs. If we denote the unit vectors In this 
coordinate system by ( 6 ] f ft?, § 3 ), then the veloc- 
ity vector can be written as 


Vi 


+ V y e 2 * V z e 3 


( 4 ) 


We now approximate the wing by an equivalent yawed 
Infinite wing at each spanwlse station. For this 
purpose, we define a new Cartesian coordinate sys- 
tem (x,y,z) constructed so that x Is normal to the 
leading edge of the wing and y Is normal to the 
airfoil chord c n . If we let a* represent the local 
incidence angle of the chord c with respect to the 
x-axis, as shown In Fig. 2, with a and b denoting 
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the leading edge of the chord line, then It can be 
shown that for an airfoil section defined normal 
to the leading edge (x/c n , y/c n ) the relationship 
between the (x,y) and (x.y) coordinates can be 
expressed by 

X a c n [(X/C n ) cosa n > (y/c n )s1na n ] 
c * c + c cosx 


i 

c 


— [— - 7- Slna * jf- COSa ] 

c c c n c _ n - 1 


(5b) 


= v • 1 * V- COSa n COSX - V- Sln<* n - V- COSa n SlnX 

(8a) 

V y = V • j = V- s1na n cos\ + V- cos<x n - V- s1na n slnx 

(8b) 

V z - V . k = V- slnx + V- cosX (8c) 

Similarly, for a freestream velocity given by 

V - V [e, cosa + e 0 slna] (9) 

a> ® l c 

the components of the freestream velocity In the 
(x,y,z) directions are: 

(V ) - V • 1 = V [cosa cosa cosx - slnaslna] 

' x ® « «> n n 

noa) 

(V ) = V • J = V [coso slna cosx + cosa slna] 

' y 7 od m J oo 1 - n n 

(10b) 

(VJ = V • k = V [cosa slnx] (10c) 

2 tt GO CD 

Since our model Is for a yawed Infinite wing, then 
V z given by Eq. (8c) and (V^*, by Eq. (10c) must be 
equal. This leads to the relation 

V- * V tanx cosa - V- tanx (11) 

z ® x 


Here a corresponds to the local Incidence angle 
of the Ohord c n with respect to the x n -ax1s, which 
represents the coordinate normal to the leading 
edge defined In the (x,z) plane. a n Is related to 
and X by 


tana^ 
tana n * cosx 


(6a) 


The relation between the. airfoil chord length nor- 
mal to the leading edge of the wing and the stream- 
wise airfoil chord length Is 


COSanCOSX 

c n * C0Sa n 


(6b) 


The trajectory calculations require that the veloc- 
ity components In three directions are given and, 
since the jjanel method computes the velocities for 
the (x,y,z) coordinate system. It Is necessary to 
determine them In the (x,y,z) coordinate system so 
that Ice accretion can be computed on an equivalent 
two-dimensional model. If we denote the un1t A vec- 
tors of the coordinate system (x,y,z) by (i,j,k), 
then the relationship between the unit vectors In 
the two coordinate systems are 

1 « e 1 cosa^cosx - e 2 s1na n - e 3 cosa n s1nX 


J = e 1 s1na n cosx + e 2 cosa n 


e 3 s1n ° n s1nx - (7b) 


slnx + e 3 cosx 


(7c) 


Each velocity component, V x , V y , V,, In the (x,y,z) 
coordinate system can be obtained by taking a dot 
product of the velocity vector V with Its respec- 
tive unit vector, that Is, 


which allows the expressions given by Eq. (8) to 
be expressed In the following forms 

2 

COSa COSa COSa $1n X 


V » V- 7^ 

X X COSX 


- V- Slna - V 
y n « 


cosx 


(12a) 


slna slna COSa s ^ n 

V y * v x cosx + V y C0Sa n " V ® cosx 


V, * V cosa slnx 
z ® 


(12b) 
(12 c) 


with total velocity given by, 

Vj » (V 2 ♦ V 2 + V 2 ) = V? t V? 

, (V^cosa - V-) 2 tan 2 X (13) 


The trajectory calculations for an airfoil 
section In a plane normal to the leading edge of 
the wing are similar to those for a two-dimensional 
flow. Instead of solving two equations of motion 
In the x and y directions, we now need to consider 
three equa tlons In the x , y and z d 1 rec 1 1 ons even 
though the particle paths are calculated only Tn a 
plane normal to the leading edge of the wing 
according to the model discussed previously. To 
elaborate on this point further, we use the tra- 
jectory analysis discussed In Ref. 3 and apply 
appropriate modifications to Incorporate our equiv- 
alent Infinite swept wing approach Into the two- 
dimensional analysis. 3 

Let us consider a water particle In a flowfleld 
where the freestream velocity makes an angle p 
with the x-axIs, as shown In Fig 3. The (x.y) 
coordinate system represents the airfoil coordinate 
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Fig. 3. Water particle velocities. 

system with \L denoting the local total velocity to 
which the particle Is subjected and is computed by 
the panel method. The particle velocity Is denoted 
by Vp and the relative air velocity that the parti- 
cle experiences by 


where 


px 


9 P 


(14) 


V 


(15a) 

dx„ • 

dy^ * 


— & i 

dt 1 

♦dt*J 

(15b) 


From Eqs. (14) and (15b), the relative air veloc- 
ity with respect to the moving particle can be 
expressed as, 

9 . (V - V )1 ♦ (V - V )j (16) 
r 1 x px' ' y py /J ' 

so that the unit vector, e , In the ? direction 
can be computed by 1 


Resolving Eq. (18) Into Its components and 
using the solutions In Eq. (19), we can write 

^ x p c <J pA P 

V V v x - V * g s1np (20a) 

- “ST 1 V\ ' V ‘ • C0S “ (Zm) 


^ z p ^pA 

^2 - V V V Z - V pz ] < 2(k > 

where 

dz^ 


If we now Include the components of the particle 
and Invlscld air velocities In the k-dlrect^on In 
the definition of the relative air velocity V r with 
respect to the moving particle, we can write the 
resulting relative air velocity as 

9 r ' <*, - V’ * <V > ' V J ’ * <»I - <“> 

If we now Invoke the yawed Infinite swept wing 
assumption which requires that the spanwlse Invls- 
cld velocity component V 2 Is constant, then we 
can write 

V z = slnx cosa (23) 


With this assumption, we note that with A * 
C c jpA p /2m, Eq. (20c) Is of the form 

dV 

= AV dt = d6 (24) 

z " pz 


since 


dV 

dt 



Assuming that the right-hand side of Eq. (24) Is a 
function of t only, we can Integrate Eq. (24) and 
write the resulting expression In the form 


e 


r 




(17) 


We now apply Newton's law of motion to the particle 
In the e r -d1 rectlon. Neglecting the lift force, 
the forces acting on the particle correspond to the 
drag force 4 F a , A and^gravltatlonal force, Fg, so 
that with r p = x) + yj and with m denoting the par- 
ticle mass, 

* f . (,8 > 


where 


; rV r A P 


(19a) 


V - V z - < V z - V Jz )e ' 6(t) (25) 

where denotes the particle release velocity. 
Equation*^ 25) shows that If the particle Is re- 
leased with a speed equal to the total Invlscld 
velocity, then the spanwlse velocity component of 
the particle Is equal to the Invlscld spanwlse vel- 
ocity which Is constant. As a result, the relative 
air velocity for a yawed Infinite swept wing, Eq. 
(22), Is the same as that for a two-dimensional 
flow, Eq. (20). This means that the expressions 
used In trajectory calculations for two-dimensional 
flows can also be used for finite swept wings when 
the flowfleld Is computed by a three-dimensional 
panel method and the spanwlse velocity component Is 
assumed to be fixed by a yawed wing approximation. 


f * mg( 1 slnp - J cosp) (19b) 


4.0 Three-Dimensional Interactive Boundarv-Laver 
Method 


Here Is the drag coefficient, p the density 
of air, Ap the cross-sectional area perpendicular 
to the flow, and g the acceleration of gravity. 
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As In two-dimensional flows, the Interactive 
method for three-dimensional flows Is based on the 
solutions of the Invlscld and boundary-layer 


equations. An Interface program, Illustrated by 
Fig. 4, Is placed between the Invlscld and three- 
dimensional Inverse boundary-layer methods to pro- 
the geometry and Invlscld velocity data for 
Input to the boundary-layer program. The basic 
Input to this program Is (1) the definition of the 
wing configuration that Is used by a geometry sub- 
routine to construct a nonorthogonal coordinate 
system and (2) the associated geometrical param- 
eters, which consist of the geodesic curvatures and 
metric coefficients needed In the boundary-layer 
calculations. Some of the generated data Is used 
later In a velocity subroutine to determine the 
Invlscld velocity components at the boundary-layer 
grid points and to transform the Invlscld velocity 
components on the surface, calculated In a 
Cartesian coordinate system, Into the boundary- 
layer coordinate system. This operation consists 
of calculating dot products of velocity vectors as 
well as chordwlse and spanwlse Interpolation. 
Further velocity and geometry data processing Is 
carried out In a subroutine that separates the 
generated Information Into upper and lower surfaces 
of the wing for boundary-layer calculations. 



Fig. 4. The Interactive boundary-layer method. 

The above procedure Is appropriate to wings 
without Ice and has been used to compute transonic 
flows on wing/body configurations 10 where, since 
the wing leading edge was free of Ice, there was 
no difficulty In generating solutions near the 
attachment line by constructing the nonorthogonal 
coordinate system and computing the geometrical 
properties of the wing. For a wing with Ice, gen- 
eration of the boundary-layer solutions near the 
leading edge can pose problems since the geodesic 
curvatures and metric coefficients must be deter- 
mined for an Irregular surface. In addition, the 
formulation of the Interactive boundary-layer 
method developed for Iced airfoils must take 
account of the three-dimensional nature of the 
flow. Thus, It Is necessary to make changes In 
the strategy for solving the three-dimensional 
boundary-layer equations for an Iced wing. These 
are considered below. 

4.1 Boundary-layer Equations 

The three-dimensional boundary-layer equations 
for a nonorthogonal coordinate system are given In 
several references (for example, Ref. 26). With 
Reynolds stresses^ modeled by the eddy-viscosity 
concept, they can be written as. 


K 0 w 2 cote + K,u 2 cosece 
h 1 ax h ? 32 ay 2 1 

2 

cote cosece ap cosec e ap 
+ K 21 UW * ^ ax ' « 


+ 


a /K aw. 
v dy 3y* 


(28) 


Here x denotes the coordinate along the lines 
formed by the Intersection of the wing surface and 
planes representing constant percent semispan; z Is 
the coordinate along the constant percent chord- 
lines that generate the wing surface, with chord 
defined as the maximum length line between leading 
edge and trailing edge. The third coordinate y 
denotes the direction normal to the wing surface, 
and the parameter h denotes the metric coeffic- 
ients, with 0 the angle between the coordinate 
lines x = const and z = const. For an orthogonal 
system, 0 = */2. The parameters K-j and K 2 are 
known as the geodesic curvatures of the curves z ^ 
const and x - const, respectively. Equations (26) 
to (28) are subject to the following boundary con- 
ditions 


y a 0, u * 0 , v * 0, w - 0 (29a) 

y a 6, u 3 u e (x,z), w = w e (x,z) (29b) 

The solution of the above equations also 
requires Initial conditions on two Intersecting 
planes; those In the (y,z) plane at a specified 
chordwlse location are determined from the solu- 
tions of the equations discussed In Subsection 4.3. 
Those on the ( x, y) plane, at a specified spanwlse 
location z = z 0 , with z 0 corresponding to, say, the 
root location, are determined from the solutions 
of the quasl-three-dlmenslonal form of Eqs. (26) 
to (28) with all derivatives with respect to z 
neglected, that Is, 

(uh 2 s1ne) «■ (vh^ sine) = 0 (30) 

, v I 51 - K.u 2 cote ♦ K_w 2 cosece ♦ K.uw 
ax ay 1 2 i 

_ cosecfe 3£ , 3_ (b M) (31) 

P h 1 3x 3y ' 3y' 

h~ 3x * V ay ” K 2 w2<:0te * K -| u2cosece * K 21 UW 
, , tote SQS ece g , v {_ ( b ft ( 32 ) 

subject to the same boundary conditions given In 
Eq. (29). 


4.2 Interaction Law 


(uh 2 sine) «■ §y (wh, sine) + ( vh ] h 2 sine) 


(26) 


U_ lit , w_ IK * v M . k n 2 cote ♦ K,w z cosece 
h 1 ax h 2 az ay l 2 

„ cosec 2 e 3fi . cote cosece 3fi 

♦ k 12 uw ‘ “ ph 1 ax ph 2 3z 

3 / h 

* v a7 (b ay* 


( 27) 


To account for possible flow separation, as In 
two-dimensional flows, we use the Interaction law 
of Veldman 27 where, for airfoil flows, the edge 
velocity Is expressed as the sum of an Invlscld 
velocity u{2(x) and perturbation velocity 4u e (x) due 
to viscous effects, that Is, 

u e (x) = ug(x) ♦ 6u e (x) (33) 

The perturbation velocity Is given by the Hilbert 
Integral 
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“V*) ■ l'* b h <V*> 

x a 


In the Interaction region (x a , x b ). 


(34) 


To extend this Inverse formulation to three- 
dimensional flows, It Is necessary that the two- 
dimensional Interaction formula given by Eq. (34) 
be either modified to account for the Interaction 
In the x- and z-dlrectlons or be replaced by 
another formulation which provides a relationship 
between displacement surface and external velocity. 
Here we use the former approach, as described In 
Ref. 18, and first generate an Initial displacement 
surface by solving the quasl-three-dlmenslonal 
boundary-layer equations subject to the boundary 
conditions given by Eqs. (29) and (33) with the 
external velocity distribution u§(x) obtained from 
the panel method. The second step Involves Inter- 
action between the Invlscld flow equations! and the 
quasl-three-dlmens lonal flow equations. As In two- 
dimensional flows, the solutions of the boundary- 
layer equations are used to compute distributions 
of blowing velocity on the surface, and these allow 
the Invlscld flow solutions to be updated. In step 
three, after the calculation of the Initial condi- 
tions In the (y,z) and (x,y) planes, the fully 
three-dimensional boundary-layer equations are 
solved with the external velocity components 
resulting from step two. As before, the spanwlse 
velocity component Is assumed to correspond to Its 
Invlscld value. The viscous flow solutions are 
obtained by marching In the spanwlse direction at 
each advancing chordwlse location. This represents 
the first phase In an Interactive loop that In- 
volves the fully three-dimensional boundary-layer 
equations. In the subsequent phases, as before, 
the blowing velocity distribution Is used to obtain 
Improved Invlscld flow solutions, so the fully 
three-dimensional boundary-layer equations are 
solved for Iced wings as for clean wings In tran- 
sonic flow.™ 


vh 1 h 2 s1ne - - (|* + |f) (36) 

In addition, dimensionless parameters f and g are 
defined by 

1 > = (u e vs) 1/2 h 2 s1nef(x,z,n) 

(37) 

l/? % 

* = (u g vs) — h.,s1neg(x,z,n) 


(bf ") ‘ + ef* + m 2 (f ') 2 ♦ mjf'g' ♦ m 0 (g ' ) 2 + ra^ 


» m in f ♦ m 7 g‘ 


"10 


ax 


az 


(38) 


(bg')' ♦ eg" ♦ n^f'g' ♦ n» 3 (g ' ) 2 t m 9 (f') 2 + i" 12 


SC- ♦ m,g‘ ?a! 


1 io f * ax - * IB 7 9 ’ az 


e- . roi f' ♦ m 6 g' ♦ m 10 fr ♦ ”7 If 1 

The coefficients m-, to m 12 are defined by 
m 1 = (zu^s h 2 s1ne) -1 ^ |j (zu^s h 2 slne) 


(39) 

(40) 


au u_ 

s e . _o 


“2 u e h 1 ax 


- * sKpote, m 3 • sK 2 jp cote. 


n 4 “ _sK 2T m 5 “ ■ u 2 h 2 az SK 12 u e ' 
e 


(41) 


_ Ug 

\ - (/u^s h^lne)’ — (/u e s h^lne — ) 


The viscous effects In the spanwlse component 
w e are assumed to be second order, although their 
neglect Is contrary to the Irrotatlonallty condi- 
tion. However, trial calculations Involving vari- 
ations of both velocity conditions showed that 
errors were smaller than those associated with the 
convergence of the solutions. 


su Q u Q 2 

m 7 * ra 8 * - sK 2 ( U; ) csc6 ' 


m Q - -sK 1 jp csce, m 


10 * h 


4.3 Transformed Equations 


As In two-dimensional flows, we express the 
boundary-layer equations In transformed variables 
for computational purposes. At first, when the 
equations are solved for a prescribed external 
velocity distribution (standard problem), we use 
the Falkner-Skan transformation and a modified 
version of this transformation for the Inverse 
mode. In the former case, the Independent vari- 
ables are defined by 


u 1/2 

x * x, z - z, n » {“) y, 


J X h dx ' 

0 1 (35) 


For the dependent variables u, v and w, we Intro- 
duce a two-component vector potential such that 

ut^slne » JJ, whjSlno * |j. 


m ll 


w w 2 

[m 2 ♦ » 5 f ♦ m fl (f) ]. 
0 0 



♦ ra 13 


"la 


|Vo % 

U^ dX 
0 


w ft aw a 
m _e e 

m 7 u 2 az 

0 


In terms of transformed variables, boundary condi- 
tions given by Eq. (29) become 

n * 0: f-g-f'-g' -0 

(42) 
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The form of the transformed quasl-three-dlmenslonal 
equations (Eqs. (30) - (32)) Is Identical to the 
form of Eqs . (38) to (40), except that 


af' ag' _ 3u e _ aw e _ Q 

az az az az “ 


(43a) 


and with w e corresponding to Its Invlscld value, 
are 


n = 0: f , g , f' , g 1 , 0 

t. = n e : f = u e . g' - w e 


(49a) 

(49b) 


and 1115 and 1% are now defined by: 

m fi » -Ki,s |r • "V = 0 <*3b) 

5 \c U g 0 

To solve the equations In the Inverse mode, we 
define Independent variables by 

u 1/2 x 

x = x, z = z, Y = (-7) y, s r j h,dx (44) 

vs 0 . 

and relate the vector potentials and ♦ to f and g 
by 1 

* = (u Q vs) 1/2 h 2 s1nef(x,z,n) 

(45) 

* » (u Q vs) 1/2 h^lnegtx^.n) 


The quasl-three-dlmenslonal form of the equa- 
tions, which are subject to the boundary conditions 
given by Eq. (49), are obtained from the above 
equations by setting 


. , . , au aw 

af 1 e e 

az * az * az az 


and m 6 - 0 (50) 


To generate the Initial conditions near the 
leading edge of the Iced wing, we use quasl-three- 
dlmenslonal boundary-layer equations expressed In 
the Inverse mode given by 

(bf)' ♦ ef’ + m 2 [(f') z - (u e ) 2 ] ♦ m 5 (f'g' - u e w e ) 


♦ m 8 t(g') 2 - (w e ) 2 ] 


«W f ' 


afj 

ax 


- au e 
- u e i 

(51) 


and with a prlme^now denoting differentiation with 
respect to Y and Up and w e denoting edge velocity 
components normalized by reference velocity u 0 , 
Eqs. (26) to (28), with e‘ defined by Eq. (40), are 
written as 

(bf-) ’ ♦ ef" + m 2 [(f') 2 - (u e ) 2 ] ♦ m 5 [f'g' - u e w e ] 


♦ m 8 I ( 9 ’ ) 2 - (w e ) 2 ] = n. 10 ( f * fj- - u e 5 f> 


au 

. 211 . n A 


af' - 3 °e 
* V 9 ' az " w e az * 

(46) 


(bg-)' ♦ eg* ♦ m [(g 1 ) 2 - (wj 2 ) ♦ m^f'g' - 0^3 


dw 

♦ m 9 [(f) 2 - (u e ) 2 ] = m 10 (f' {£■ - u. -^r) 


e ax 


, , aai - A 

+ m 7< 9 az " w e az _) 
(47) 

The coefficients m-, to m 1Q are now given by 

5 1/2 a 1 /p 

"l * h 1 h 2 sine ax <s h 2 s,ne) ' m 2 * sK^oto, 
m 3 = $K 2 cot8, m^ 3 -sK 2 ^» m^ - 

*. - f; <’ ,/ V<"'» <«> 


m 7 * s/h 2 , m 8 = -sK 2 csce, 


m^ * -sK-jCSC©, m-jQ * s/h-j 

The transformed boundary conditions for the system 
of Eqs. (46) and (47), with u e given by Eq. (33) 


(bg")' * eg* + m 3 ((g') 2 - (w,) 2 ] ♦ m 4 <f'g' - u^) 

3w 

♦ ra 9 [(f) Z - (0 e ) 2 l * m 10 (f g 1 - u e j-S) 

(52) 

e‘ = m,f 1 ♦ m. n (53) 


The above equations can be further simplified 
If we assume that two adjacent defining sections 
of a wing are connected by straight line develop- 
ment, as commonly used In the wing design. This 
feature simplifies the problem of shaping the metal 
for a wing surface. As a consequence, we can 
neglect the geodesic curvature of x = constant 
lines, namely K?, and thus set m3 = mg * 0. From 

the definitions of 1114 and m$, It can be seen that 
as a result of the above assumption, these two 
terms are also small and can be neglected. We 
further assume that the local variations in cross 
sections In the spanwlse direction are small. 
Examination of the terms m] , m2 and mg for a typi- 
cal wing shows that m2 reaches a value less than 
0.1 very close to the leading edge (x/c < 0.01) 
and mg reaches a maximum value of 0.2. However, 
their magnitudes rapidly decrease with Increasing 
x and reach a very small value at x/c < 0.1. This 
behavior allows us to neglect m2 and mg In the 
equations and set m^ * 1/2. 

4.4 Solution Procedure 

A detailed description of the solution proced- 
ure will be reported later separately. Briefly, 
the boundary-layer equations expressed In terms of 
transformed variables are solved with Keller's two- 
point finite-difference method 2 ^ (box scheme) 
with boundary conditions expressed In Inverse form 
with the Interaction law described In Subsection 
4.2. Depending on the complexity of the flowfleld, 
two forms of the box scheme are employed. In 
regions where all velocity components are positive, 
the regular box scheme Is used. In regions of 
either a negative spanwlse velocity component or 
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negative streamwlse velocity component, the zig-zag 
box scheme described In Ref. 26 Is used. 


5.0 Comparison of Calculated and Measured 
Ice Shapes on a Swept Wing 

The approach described In Section 3 was applied 
to the calculation of Ice shapes that have been 
observed on a MS-31 7 swept wing obtained under 
Infinite swept wing conditions. These data were 
also used for a similar purpose by Potapczuk and 
Bldwell^, as discussed In Section 3. 

In general, the Infinite yawed wing conditions 
apply to the mld-semlspan section of wings with an 
aspect ratio greater than about five. This approx- 
imation becomes progressively less accurate as the 
tip or the root of the wing Is approached but In 
most Instances It can still provide reasonable 
answers. A point to remember about the use of this 
approximation with finite aspect ratio wings Is 
that although the flow may have the desired char- 
acteristics, Its lift Is always less than the lift 
of a wing with Infinite aspect ratio. This may 
lead to problems In comparing calculations with 
experimental data unless the aspect ratio or the 
pressure distribution Is also given. If the pres- 
sure distribution Is not available, the given angle 
of attack may not properly represent the experi- 
mental conditions. Similar problems may also 
arise In simulating wind-tunnel conditions by 
calculating the corrected Incidence and lift coef- 
ficient In free air, because the trajectories In 
the two cases may be far from Identical. One solu- 
tion to the wind-tunnel problem, which may be the 
only acceptable solution for a swept wing spanning 
the tunnel. Is to calculate the flowfleld about the 
wing In the presence of the tunnel walls. This Is 
within the capability of the Invlscld method dis- 
cussed In Section 2, and studies are underway using 
this capability to Investigate the effect of wind 
tunnel walls on the pressure distribution. 

The comparisons presented here are for the test 
conditions given In Table 1. Additional studies 
for other test conditions are In progress and will 
be reported later. The calculated Ice shapes In 
Figures 5 and 6 were obtained for the untapered 
wing with a MS-317 airfoil section defined stream- 
wise with a sweep angle of 30° and an aspect ratio 
of six. All trajectory and Ice accretion calcula- 
tions were carried out with Invlscld flow computed 
on the mld-semlspan section where the spanwlse 
pressure gradient was negligible. All calculations 
were performed for one time step to save computer 
time, which Is approximately 7 minutes per run on 
the Cray computer. The Increase In time, In com- 
parison with the two-dimensional case, Is primarily 
due to the trajectory calculations where, despite 



Fig. 5. Comparison of calculated (solid lines) and 
measured (dashed lines) Ice shapes. Rime Ice: (a) 
Run 8, (b) Run 11, (c) Run 9, (d) Run 10. 


Table 

1. Test 

Conditions 

for MS-317 

Ice Accretion 


Experiment of Ref. 

8, V* * 150 

mph, d * 20 ym 



LUC 

=*1.03 gm~ 3 



T 

t 

a 


Run 

e?) 

(sec) 

(deg) 

(k s /c)i 

1 

15 

1164 

2.0 

0.00192 

6 

15 

1164 

8.0 

0.00192 

7 

15 

390 

2.0 

0.00192 

8 

0 

390 

2.0 

0.00127 

9 

0 

390 

8.0 

0.00127 

10 

0 

1164 

8.0 

0.00127 

11 

0 

1164 

2.0 

0.00127 


the yawed Infinite wing approximations, the compu- 
tation of the off-body velocities Involves repeated 
large matrix multiplications In which all wing 
panels are represented. 

Figure 5a shows a comparison of measured and 
calculated Ice shapes for Run 8 which corresponds 
to T. » 0*, a . 2', t * 390 sec. As can be seen, 
the agreement between measured and calculated 
results Is remarkably good. The calculated results 
for a calculation time of 1164 sec and for T„ = 0*F 
and a » 2* (Run 11) are shown In Fig. 5b and Indi- 
cate reasonable agreement with measurements despite 
the one time step used In the calculations. It Is 
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6.0 Performance Degradation of an Iced 
Tapered Wing 



Fig. 6. Comparison of calculated (solid lines) and 
measured (dashed lines) Ice shapes. Mixed Ice: (a) 
Run 7, (b) Run 1, (c) Run 6, 

expected that the Ice growth will have some effect 
on the velocity field and on the calculated droplet 
Impingement. A comparison of predicted and mea- 
sured Ice shapes obtained for To, = 0°F at a = 8° 
for t = 390 and 1164 sec. (Runs 9 and 10) are shown 
In Figs. 5c and 5d, respectively. The agreement 
Is again reasonable, keeping In mind that only one 
time step was used In the calculations. 

The next set of studies was conducted for a 
slightly higher freestream temperature of T* = 
15°F, representing an Icing condition for which a 
mixed Ice growth was observed. Run 7 In Fig. 6a 
for a = 2° and t * 390 sec. Indicates good agree- 
ment between experiment and theory, except for some 
deviation on the upper surface. The results In 
Fig. 6b at the large time step of t » 1164 sec. 
(Run 1) are more or less In agreement In predicting 
the amount of Ice accumulated, but they differ In 
predicting Its shape. It Is known from two- 
dimensional calculations that a large number of 
relatively short time steps are needed to predict 
horn-shaped Ice for glaze Ice. Since the mixed 
Ice formation tends toward glaze ice shapes for 
large times. It Is not surprising that one time 
step calculation Is not sufficient to predict the 
actual growth of the Ice shapes. Similar comments 
apply to Fig. 6c, where comparisons are for a large 
time step of t * 1164 sec. (Run 6), but at o » 8°. 


The Interactive boundary-layer method of Sec- 
tion 4 was used to study the performance degrada- 
tion of an Iced wing having MS-317 airfoil stream- 
wise sections, an aspect ratio of 3.43, and a taper 
ratio of 0.4. Icing conditions were chosen to 
correspond to those In Runs 8 and 11, shown In 
Table 1. The pressure distribution on the wing was 
computed at four locations defined by the midsec- 
tion of each wing-section with a hundred panels on 
each defining airfoil section. The Ice shapes 
corresponding to this pressure distribution were 
computed with the method of Section 3 In the middle 
of each wing section and were used to distribute 
Ice along the leading edge of the tapered wing. 
The computed Ice shapes for a ^ 2° were then 
assumed to be the same for all angles of attack on 
the wing In the Investigation of the performance 
degradation of the wing due to Ice shapes corres- 
ponding to the atmospheric conditions given In Runs 
8 and 11. At a specified angle of attack, with the 
defined Ice shapes along the leading-edge of the 
wing, calculations were performed with the method 
of Section 4; that Is, Invlscld flow calculations 
performed for an Iced wing were followed by the 
Inverse three-dimensional boundary-layer calcula- 
tions to determine the blowing velocity distribu- 
tion to be used In the Incorporation of viscous 
effects Into the Invlscld method. The Invlscld 
flow solutions made use of four lifting strips, and 
the viscous flow calculations Included boundary- 
layer calculations on the wing and In the wake, the 
latter requiring velocities at off-body points In 
the potential field. This Interactive and Iter- 
ative procedure was repeated until the solutions 
converged. The lift coefficients were then calcu- 
lated from the Invlscld method for each Individual 
strip and Included the contribution of Ice protrud- 
ing beyond the wing contour and the drag coeffic- 
ients from the boundary-layer calculations. 

Figure 7 shows the variation of the calculated 
lift coefficients as a function of angle of attack. 
Since the primary purpose of the calculations was 
to demonstrate the Increase In drag due to Ice on 
a tapered wing, the angle of attack range was not 
extended to stall, which would occur at relatively 
high angles of attack for low aspect ratio wings. 
The higher lift coefficient than for the clean wing 
shown for the two Iced wings Is due to the normal- 
ization with the wing area of the clean wing In 
both cases. The conclusion from this figure Is 
that lift Is not affected by the rime Ice accretion 



Fig. 7. Effect of leading-edge 390 and 1164 second 
rime Ice on the lift coefficient of a tapered wing 
for R = 4.6 x 10 6 based on root chord. 


for the angle of attack range considered here 
because the Ice shapes along the leading edge of 
the wing for runs 8 and 11 do not cause premature 
flow separation on the wing. 

The calculated drag coefficients shown In Fig. 
8 represent the profile drag of the wing only and 
do not represent the total drag, since that 
requires the contribution of the Induced drag. The 
profile drag was calculated sectlonwlse from the 
Squire-Young formula based on the resultant veloc- 
ity at the trailing edge. Comparable results were 
also obtained from the momentum deficiency In the 
far wake. Here we see considerable differences 
between the clean wing and the two Iced wings 
because the Reynolds number Is relatively low (Re 
= 4.6 x 10^ for the root chord) and there are 
large regions of laminar flow on the clean wing. 
The principal contributor to the drag' Increase for 
the Iced wing Is the shift In transition to near 
the leading edge due to roughness of the Ited sur- 
face. The contribution of the surface roughness 
Itself to the drag Is very small for Run 8 because 
the extent of Ice Is small and Its shape emulates 
an airfoil leading edge. The additional drag 
Increase for Run 11 results from the surface rough- 
ness spread over a large wetted area Increment. 
The main conclusion that can be drawn from these 
comparisons Is that drag Increments obtained be- 
tween clean and Iced airfoils In wind tunnels 
depend on transition locations on the clean wing. 
If the Run 8 case represents a wing with transition 
fixed at the leading edge and the clean wing case 
Is transition free, the observed drag Increments 
from the Run 11 case are quite different from each 
other. As a corollary, drag Increments obtained 
from wind-tunnel tests may be meaningless without 
fixing transition or knowing where transition 
occurs during the tests. 



Fig. 8. Effect of leading edge 390 and 1164-second 
rime Ice on the profile drag coefficient of a 
tapered wing for R c = 4.6 x 10^ based on root 
chord . 

7.0 Concluding Remarks 

A method for predicting Ice accretion along the 
leading edge of a wing Is described and evaluated 
by comparing calculated results with experimental 
data for several atmospheric Icing conditions ob- 
tained on a MS-317 wing with Infinite swept-wlng 
conditions. Overall, calculated Ice shapes are In 
good agreement with measurements. 

A method Is also described for computing the 
aerodynamic properties of clean and Iced wings. 
This method Is based on the extension of the 


Interactive boundary-layer method previously devel- 
oped for clean and Iced airfoils. It Is applied 
to a tapered with wing leading edge rime Ice accre- 
tion of 390 and 1164 seconds to study the wing 
performance degradation for a range of angles of 
attack less than stall at a root-chord Reynolds 
number of 4.6 x 10°. Calculated lift coeffic- 

ients Indicate that the lift Is not affected by the 
rime Ice accretion because the leading-edge Ice 
shapes do not cause premature flow separation on 
the wing for the atmospheric Icing conditions and 
angle of attack range considered In the study. 
Calculated profile drag coefficients for clean and 
iced wings, on the other hand, show the expected 
differences; that Is, the Increase In drag due to 
leading-edge Ice. The main reason for the drag 
Increase Is the movement of the transition loca- 
tions to the wing leading edge due to the roughness 
of the Iced surface. The contribution of the sur- 
face roughness to the drag Is relatively small for 
the small ice accumulated for 390 seconds because 
the extent of Ice Is small. For an Ice shape 
corresponding to 1164 seconds, however, the drag 
Increase Is more pronounced because. In this case, 
the surface roughness spreads over a larger wetted 
area Increment. 

Both methods for predicting Ice shapes along 
the wing leading edge and for computing Its aero- 
dynamic performance degradation due to Icing are 
very encouraging. The methods are general and can 
be applied to other parts of an airplane other than 
a wing. Including the engine Intakes. Studies are 
currently underway to apply these methods to deter- 
mine the performance degradation of a typical twin- 
engine commuter-type aircraft In measured natural 
Icing conditions. 
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